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ABSTRACT 


Optimal  control  algorithms  that  use  an  adaptive  non-recursi ve 
digital  filter  model  for  on-line  closed-loop  blood  pressure  regulation 
have  been  developed.  This  Automatic  Therapeutic  Control  System  was 
designed  specifically  for  the  regulation  of  physiological  systems,  but 
the  design  assumptions  are  such  that  it  should  prove  very  useful  for 
a  much  broader  class  of  control  problems. 

An  Adaptive  Model  Control  system  has  been  developed,  the  analysis 
of  which  is  presented  in  two  parts:  (1)  the  real-time  adaptive  model 
synthesis  procedure,  and  (2)  the  optimal  forward-time  controller. 

The  Adaptive  modeling  process  is  accomplished  by  the  rapidly 
converging  O-LMS  Algorithm.  The  conditions  necessary  to  guarantee 
convergence  for  deterministic  inputs  are  presented.  Mean-Square  Error 
bounds  are  presented  for  zero  mean  and  nonzero  mean  additive-output 
noise  systems  and  also  for  the  low-order  approximation  problem. 

The  optimal  forward-time  controller  is  described.  It  not  only 
makes  efficient  use  of  the  mathematical  properties  of  the  non-recursi ve 
digital  filter  model  (that  is,  the  filter  gains  and  the  state  of  the 
filter),  but  also  meets  the  time  and  memory  constraints  of  a  mini¬ 
computer  used  on-line  for  this  control  problem.  The  future  control 
inputs  are  determined  by  a  doubly  constrained  quadratic  function  which 
is  solved  to  minimize  the  mean-square  control  error. 

The  results  of  an  experimental  run  are  included  in  which  these 
algorithms  were  used  to  regulate  the  blood  pressure  of  a  dog  that  had 
been  artificially  placed  in  a  hypotensive  state  (shock).  The  results 
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of  this  and  similar  experiments  have  been  very  successful  from  both 
an  engineering  and  medical  point  of  view,  and  if  the  necessary  arrange¬ 
ments  can  be  made  with  the  local  hospitals,  this  system  will  be  used  in 
the  near  future  as  part  of  an  intensive  care  unit. 
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1.  INTRODUCTION 


A.  AUTOMATIC  THERAPEUTIC  CONTROL  SYSTEMS 

Therapeutic  is  defined  as  that  which  is  related  to  the  branch  of 
medicine  dealing  with  the  treatment  of  disease  [l].  Thus,  an  automatic 
therapeutic  control  system  is  any  control  system  which  deals  with  the 
treatment  of  an  illness,  sickness,  interruption  or  perversion  of  function 
of  any  of  the  organs  in  an  automatic  (closed  loop)  manner.  The  automatic 
therapeutic  control  system  proposed  in  this  paper,  an  adaptive  model 
controller  (see  Fig. 7),  was  used  for  the  treatment  of  shock  (hypotension). 
For  this  problem  a  pressure-elevating  drug  was  used  to  regulate  the 
average  blood  pressure  of  a  dog  which  had  been  placed  in  a  controlled 
state  of  shock. 

Another  promising  area  of  application  is  anesthesiology.  Many  of 
the  chemicals  used  in  anesthesia,  both  in  the  liquid  and  in  the  gaseous 
form,  are  sufficiently  fast  acting  to  utilize  fully  this  control  system. 
This  system  for  quickly  controlling  the  degree  to  which  a  patient  is 
anesthetized  would  thus  be  a  great  help  to  an  anesthesiologist. 

B.  THE  CARDIOGENIC  SHOCK  PROBLEM 

It  is  estimated  that  in  the  United  States  alone  250,000  lives  are 
lost  annually  because  of  cardiogenic  shock  [2,3],  The  best  treatment 
for  this  deadly  problem  is  still  disputed,  and  it  is  claimed  that  cardio¬ 
genic  shock  still  defies  the  cardiologist's  skill  80  to  90  per  cent  of 
the  time. 

Thus,  there  is  a  need  for  new  and  improved  ways  to  deal  with  cardio¬ 
genic  shock.  One  possible  method  of  treatment  is  to  regulate  automati¬ 
cally  a  patient's  blood  pressure  by  administering  drugs  with  a  control 


system.  The  engineering  control  problems  associated  with  this  approach 
are  what  we  are  concerned  with  in  this  paper. 

The  methods  of  treatment  for  cardiogenic  shock  can  presently  be 

broken  into  two  opposing  schools  of  thought:  those  who  prefer  drug 

therapy  and  those  who  prefer  mechanical  intervention,  particularly  with 

the  intra-aortic  balloon  pump.  We  are  concerned  with  the  drug  therapy 

approach,  which  can  be  broken  into  three  schools  of  thought  [2]. 

"The  vasopressor  school,  which  is  becoming  popular  again, 
emphasizes  the  need  for  agents  which  increase  coronary  per¬ 
fusion  pressure.  But  proponents  of  vasodilators  contend 
that  the  major  need  is  to  relieve  intense  vasoconstriction. 

This  mode  of  therapy  reduces  the  work  load  of  the  heart  and 
redistributes  blood  flow  to  all  areas  of  the  body.  Still  a 
third  group  emphasizes  the  need  for  increased  myocardial 
contractility  to  restore  cardiac  output." 

Regardless  of  the  approach,  the  objective  is  to  take  therapeutic 
measures  early  enough  to  prevent  irreversible  shock  from  developing. 
Irreversible  shock  can  be  described  as  a  state  of  positive  feedback 
where  the  failure  of  one  function  causes  still  further  failure  of  another 
function  within  the  same  control  loop  [ 4 , 5 ] .  Figure  1  is  a  simplified 
block  diagram  of  some  of  the  different  types  of  feedback  that  can  lead 
to  a  progressive  state  of  shock  (irreversible  shock)[5]. 

We  again  remind  the  reader  that  this  paper  is  concerned  with  de¬ 
scription  and  analysis  of  a  control  system  which  regulates  blood  pres¬ 
sure;  we  do  not  address  the  larger  problem  of  what  therapy  is  best  for 
cardiogenic  shock. 

C.  THE  COMPUTER  LABORATORY  &  OPERATING  ROOM  FACILITIES 

The  computer  laboratory  used  in  this  research  is  located  on  the 
Stanford  Campus  in  the  Durand  Laboratory,  and  the  operating  room  is 
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Different  Types  of  Feedback  That  Can 
Lead  to  Progression  of  Shock. 


Fig.  1 . 


located  in  Palo  Alto  at  the  Palo  Alto  Medical  Research  Foundation.  A 
voice  grade  telephone  line  is  used  as  the  data  link  between  the  two 
facilities  (see  Fig.  2).  This  choice  of  data  link  will  make  it  very 
inexpensive  and  easy  to  reach  patients  in  the  intensive  care  units  at 
any  of  the  local  hospitals. 

The  primary  hardware  components  are  shown  in  block  diagram  form  in 
Fig.  3.  In  the  operating  room,  the  strain  gauge  and  drop  master*  are 
commercially  available  items,  while  the  voltage-controlled  oscillator 
(VCO)  and  drop-rate  DETECTOR  were  designed  and  built  at  Stanford.  The 
acoustic  couplers  and  the  HP-2116B  computer  system  are  commercially 
available  items,  but  the  frequency-to-digital  converter  (FDC)  and  the 
drop-rate  ENCODER  were  also  designed  and  built  at  Stanford. 

The  frequency  band  from  lKHz  to  3  KHz  is  used  for  transmission  of 
the  instantaneous  blood  pressure.  The  remaining  channel  space  of  the 
telephone  line  (500  to  1000  Hz)  is  used  to  transmit  the  encoded  drop 
rate  commands.  A  tone-burst  code  is  transmitted  each  time  a  drop  of 
drug  is  required. 

The  HP-2116B  computer  is  a  mini-computer  with  a  16,000  word  (16 
bits  per  word)  memory  and  a  1.6  microsecond  cycle  time.  In  addition  to 
the  major  components  shown  in  Fig.  3,  this  computer  system  has  a  real 
time  clock  which  is  required  for  scheduling  data  gathering,  data  process¬ 
ing  (the  forward  time  calculations  and  adaptive  algorithms)  and  data 
output . 

*A  drop  master  is  a  machine  which  releases  one  drop  of  fluid  in  response 
to  either  an  internal  timer  or  an  external  command. 
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Pictorial  Representation  of  the  Experimental  System 


BLOOD  PRESSURE 


Fig.  3.  Block  Diagram  of  the  Primary  Hardware  Components. 

Interested  readers  are  referred  to  references  [6-8]  for  a  more  com¬ 
plete  description  of  the  biomedical  engineering  problems  and  requirements 
associated  with  patient  monitoring  and  automatic  therapeutic  control 
systems . 


D.  THE  ADAPTIVE  MODEL  CONTROLLER 

When  designing  controllers  for  real  physical  systems,  we  can  cope 
with  only  those  problems  that  can  be  foreseen  and  allowed  for  in  advance. 
Thus,  the  idea  of  an  adaptive  control  system  which  can  compensate  for 
the  time-varying  and  nonlinear  aspects  of  a  system  is  very  attractive. 

The  adaptive  control  principle  in  essence  consists  of  three  things  [9]: 

1.  The  definition  of  an  optimum  condition  of  operation. 

2.  The  comparison  of  the  actual  performance  with  the  desired 
performance . 
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3,  The  adjustment  of  system  parameters  by  means  of  closed- 
loop  operation  so  as  to  drive  the  actual  performance 
toward  the  desired  performance. 

One  of  the  simplest  systems  which  employs  these  principles  is  the 
adjustable  compensator  controller  (see  Fig.  4).  The  adaptive  filter  is 
adjusted  to  give,  say,  a  minimum  mean-square  error.  The  problem  with 
this  system  is  that  the  desired  response  d^.  which  is  needed  to  adapt 
the  filter  is  unknown.  Note  that  if  d.  were  known,  then  the  adaptive 
filter  would  be  unnecessary. 

A  more  workable  adaptive  controller  is  the  inverse-model  controller 
(see  Fig.  5).  Under  the  appropriate  conditions  this  system  is  very  use¬ 
ful  [10].  Its  major  shortcomings  are  due  to  the  fact  that  we  must  have 

-N  -1 

an  estimate  for  the  delay  Z  and  the  inverse  model  G  must  be  well 
defined. 

Another  approach  is  the  feedback  model  controller  (see  Fig.  6) 
which  uses  a  faster-than- real-time  closed-loop  model  of  the  physical 
system  to  adapt  the  compensator.  While  this  is  certainly  attractive 
because  we  can  exercise  the  upper  loop  in  order  to  adapt  the  compensator 
without  disturbing  the  physical  system,  the  problem  lies  in  the  fact 
that  the  mean-square-error  performance  function  is  knownto  be  irregular, 
non-parabolic,  and  containing  relative  optima  [ll]. 

The  problems  inherent  in  the  above-mentioned  systems  are  overcome 
by  the  proposed  Adaptive  Model  Controller  of  Fig.  7.  In  this  system  the 
identification  process  is  achieved  by  the  rapidly  converging  O-LMS 
Algorithm,  and  the  adjustable  controller  is  a  faster-than- real-time 
forward-time  calculation  which  minimizes  the  squared  error  (r^-y^)^. 

In  this  system  the  control  is  optimal  in  a  squared-error  sense  if  the 
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Fig.  4.  Adjustable  Compensator  Controller. 


Fig.  6.  Feedback  Model  Controller. 


O 


jfc 

adaptive  model  has  converged.  Thus,  in  practice  the  system  is  run  open 
loop  until  the  adaptive  model  has  converged,  and  then  the  loop  is  closed 
through  the  forward-time  calculation. 

We  will  use  the  following  notation  to  distinguish  between  scalar, 
vector  and  matrix  variables : 

1.  a  variable  with  no  underlining  is  a  scalar, 

2.  a  variable  with  a  single  underline  is  a  vector,  and 

3.  a  variable  with  a  double  underline  is  a  matrix. 

For  example  (see  Fig.  8)  x.  is  scalar  input  at  time  i ,  and  x  is 

J  -  j 

the  vector  "input”  at  time  j. 

In  addition,  we  will  adopt  the  convention  that  the  weight  vector 
w^.  (see  Fig.  8)  is  of  length  N. 

*  ■  ———————— 

The  convergence  proofs  given  in  section  2  assume  that  the  physiological 
system  is  time-invariant. 
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2.  REAL-TIME  ADAPTIVE  MODEL  SYNTHESIS 


A.  WIDKOW-HOFF  ALGORITHM  (g-LMS  ALGORITHM) 

The  theory  of  least-mean-square  error  filter  design  has  received 
a  considerable  amount  of  attention  in  the  last  two  decades.  The  theory 
has  been  extended  from  the  problems  of  filtering  and  prediction  [12]  to 
include  those  of  system  identification,  process  control,  and  pattern 
recognition.  The  Widrow-Hoff  Algorithm  [  13 ]  (see  Fig.  8.) 


where 


w  •  , 

-J+l 


=  w  .  + 
-J 


e . 

J 


e  .  =  y  .  -  x  .  w  . 
J  J  ~j  ~J 


(2.1) 


was  originally  proposed  for  use  in  systems  where  little  or  no  a  priori 
statistical  information  is  available  and  where  memory  size  and  computa¬ 
tional  speed  are  limited.  The  theory  of  the  a-LMS  Algorithm  and  its 
many  variants  has  been  pursued  diligently  during  the  last  decade  [14-28], 
But  the  behavior  of  the  a-LMS  Algorithm  in  the  presence  of  highly 
correlated  or  dependent  inputs  has  not  been  analyzed.  The  previous 
methods  of  proof  basically  required  that  E[x^x_.]  be  nonsingular, 
and  that  the  sequential  input  vectors  be  uncorrelated.  The  need  to 
alleviate  these  assumptions  is  obvious  in  any  system  where  there  is  a 
real-time  identification  utilized  in  a  larger  control  process. 


The  behavior  of 
inputs  is  the  thesis 
deterministic  inputs 
modeling  of  discrete 


the  a-LMS  Algorithm  when  subjected  to  deterministic 
of  this  section.  Motivated  by  the  conditions  on 
discovered  by  Spain  [29] for  identification  and 
linear  systems,  we  are  able  to  derive  a  related 
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represent  the  physiological  system  to  be  modeled  by  a  non- recursive 
discrete  filter  G  of  length  N.  (See  Fig.  9.)  The  motivation  for 


Fig.  9.  Idealized  System  Where  the  Physiological  System 
is  Replaced  by  Its  Impulse  Response. 

this  representation  is  twofold:  First,  the  filter  gains  g  ,lal,2,...N 
are  equal  to  the  impulse  response  of  the  physiological  system  at  times 
(i-l)A  where  A  equals  the  time  delay  between  each  filter  tap.  Thus, 
from  our  experimental  data  we  choose  N  and  A  such  that  (N-l)A 
equals  the  total  time  of  the  impulse  response  and  such  that  A  is  small 
enough  to  give  a  "good"  piecewise-linear  approximation  to  the  continuous 
impulse  response.  Second  ,  this  representation  gives  us  a  one-to-one 
correspondence  between  the  elements  of  the  filter  G  and  the  elements 
of  the  weight  vector  w..  Thus,  the  modeling  error  at  time  j,  e., 
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»wrw«ww«i*»* 


can  be  written  as  follows 


e 


J 


T 

-  x  w 
-3  “J 


T  „  T 
-  x  .  G  -  x  .  w  . 

-J  ~  “J  “J 


=  xT(G  -  w  .) 
-J  “  -J 


(2.2) 


As  a  worst-case  analysis,  we  would  like  to  know  the  worst  that  can 


happen  to  IIg  -  w .  II  when  x.  is  arbitrary,  i.e.,  is  the  a-LMS  Algorithm 


~ J 


— J 


stable?  First,  we  consider  the  degenerate  case  of  x  3  0.  By  defini¬ 


tion,  if  x  30,  then  (2.1)  is  w  „  =  w . . 

-J  -J+l  -3 


Note  also  that  if  x.=  0, 
J 


then  by  (2.2)  e .  =  0.  The  general  result  is  that  without  making  any 
assumptions  about  the  input,  we  have: 


Theorem  1.  Given  a  non- recursive  filter  G  of  length  N  and  any 

input  x  then  the  next  weight  vector  obtained  by  a-LMS 

adaption  (0  <  a  <  2)  is  such  that 


w  -  G  <  w.  -G  I  (2.3) 

-j+i  —  -  -j  -! 


where  equality  holds  iff  e  =  0. 

Theorem  1  is  proved  in  Appendix  A. 

Theorem  1  states  that  no  matter  what  the  input  x^  is,  the  distance 
between  w  and  G  is  non-increasing.  Thus,  referring  to  Fig.  10,  we 

-j  “ 

can  see  that  the  tip  of  the  weight  vector  w  ^  will  always  be  contained 
inside  some  hypersphere  of  diameter  j|w^  -G|j  centered  on  the  tip  of  G. 
In  other  words,  (2.1)  is  stable  (convergent)  for  0  <  a  <  2. 
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C.  UNBIASED  CONVERGENCE 

Before  introducing  the  definition  of  an  N-dimensional  input 
sequence,  let  us  examine  some  of  the  more  obvious  properties  of  the 
algorithm  (2.1).  Take  the  case  where  a  -  1  and  the  input  sequence  is 

x[  =  [100  ...  0] 

x^  =  [0100  . . .  0] 

[0  ...  001] 

[1000  . . .  o] 


then  it  is  easy  to  show  that  given  any  finite  w  ,  we  have  w  =  G  lor 

~o 

all  j  >  N  +  1.  Furthermore,  if  a  £  1  but  0  <  a  <  2 ,  then  it  is  not 
1  i.  ni 

hard  to  show  that  w  =  G. 

J  -»  “  -J 

There  are  many  such  examples  that  demonstrate  the  unbiased  conver¬ 


gence  of  (2.1),  but  what  can  be  said  about  an  arbitrary  input  sequence 


So  that  we  may  deal  with  such  a  question  in  a  rigorous  manner, 
x.  be  the  present  input  vector,  an  arbitrary  vector  in  an 


Euclidean  space 


and  let  S  be  a  subspace  containing  the  previous 


N  -  1  vectors,  x  .  , ,x  .  . x .  „  ,  .  Then  the  present  vector  x 

— j-1  — j-2  — j-N+1  ~J 

can  be  uniquely  represented  in  the  form  [34] 


whe  re 


x  .  =  *i .  +  x  . 
-J  -*J  -.1 


n .  j.  S  and  x  .  >.  S 
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We  can  now  make  the  following  formal  definition  of  an  N-dimensional  input 
sequence : 


Definition. 


'V 


00 

O 


An  N-dimensional  input  sequence  is  any  sequence  of  vectors 
where  a  sequence  member  is 


x 

-J 


A 


x  . 

j-N+1 


such  that  (1)  the  sequence  of  vectors 
independent  element  vectors  for  all  j 
0  <  6  <  1  such  that 


*  has  N  linearly 
>  2N,  and  (2)  there  exists  a 


~  >  6  for  all  j  >  N 


The  restrictions  imposed  by  this  definition  will  be  relaxed  in  a  subse¬ 
quent  development. 

Consider  condition  (1)  with  N  -  2,  We  have  a  matrix,  B^,  say, 


where 


-i  ~ 


- 

XT 

— 

x  .  ,  X  . 

-j-1 

J-l  J-2 

_  —  — 

T 

X. 

x  .  x.  . 

3  J"1 

“J 

. 

-  J 

At  time  i,  the  scalar  x  is  the  only  new  input  value,  and  if  det B  ^ 0 

3 

then  we  must  have  one  of  two  cases:  if  x.  „=0  then  x  is  arbitrary, 

3~z  3 


else 


3  Xj-2 
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Condition  (2)  excludes  a  set  of  values  at  time  j  which  would  produce 

N 

vectors  in  £  whose  boundary  defines  a  circular  hypercone  centered 

at  the  origin  with  its  axis  perpendicular  to  S.  To  see  this  note  that 


,  ~  2  ,2 
x  .  +  n . 

|  -J  -\J  I 


>  6 


and  thus 


INI2  2  A  fell2 


A  constant  input  vector  is  firbidden  if  we  are  to  satisfy  the 

definition  of  an  N-dimensional  input  sequence.  It  is  obvious  that  we 

wish  to  exclude  such  a  sequence,  for  one  cannot  measure  system  dynamics 

without  perturbing  the  system.  The  fact  that  B.  must  have  rank  N 

corresponds  to  the  fact  that  there  are  N  states  in  the  filter  and  all 

states  must  be  excited  to  be  identified.  Equivalently,  the  filter  spans 
N 

E  therefore  must  have  N  linearly  independent  rows  if  it  is  to 

„N 

span  E  . 

Theorem  2.  Given  a  non- recursive  filter  G  of  length  N,  and  an  N- 

CO 

dimensional  input  sequence  (Xj?0>  the  weight  vector  w^  obtained 
by  a-LMS  adaption  (0  <  a  <  2)  converges  to  G  as  j  ->  ». 

Theorem  2  is  proved  in  Appendix  A.  The  method  of  proof  is  illustrated 
in  Fig.  10  for  N  =  2.  The  input  vectors  are  operated  on  in  the 

convergence  proof  in  groups  of  N  by  a  transformation  which  yields 


■n  .  =  X  ,  -  X  . 

JJ  ~3  ~3 


such  that  jj  is  orthogonal  to  all  N  -  1  previous  x^'s  in  the  group. 
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This  is  a  Gram-Schmidt  transformation  which  restarts  every  time  N 

vectors  have  been  processed.  Thus,  the  quantity  jlw.  -  G jj  at  any  time 

I!  0  “ I! 

j  is  bounded  above  by  the  convergence  due  to  the  Hj's‘  (This  follows 

from  Theorem  1  and  the  orthogonal  nature  of  the  tj  's.)  The  Hj's  are 

the  "new  information"  and  are  all  non-zero  due  to  the  fact  that  by  defi~ 

nition  B  ,  has  rank  N.  Note  that  the  closer  ti  is  to  x,  (the 
=j  -J 

closer  &  is  to  1)  the  faster  w^.  will  converge  to  G.  For  the 
example  considered  in  Section  2.C,  nj  =  x.  and  6=1  because  the 

-J  -J 

input  sequence  was  itself  orthogonal.  The  minimum  rate  of  convergence 
is  determined  by  the  size  of  5. 

It  is  of  interest  to  note  that  the  assumption  of  an  N-dimensional 
input  sequence  (x.)”  is  more  restricted  than  is  necessary  to  guarantee 
unbiased  convergence.  Assume  that  we  do  have  M-  N  (where  M  >  N) 


r 

dependent  vectors  in  the  sequence  { ^  J q . 


By  Theorem  1  we  know  that 


|W  -  G II  is  non-increasing.  Also,  when  a  dependent  vector  x  is  used 
to  adapt  the  weight  vector  vr  ,  the  upper  bound  defined  by  the  tj’s 
(see  Fig.  10)  does  not  change  because  jrj  is  equal  to  zero.  Thus,  the 
upper  bound  does  not  change  for  M  -  N  adaptionsbut  decreases  normally 
for  the  remaining  N  adaptions.  Therefore,  the  first  part  of  the  defini¬ 
tion  of  an  N-dimensional  input  sequence  can  be  modified  to  read  ..."such 

r  ,  j-1 

that  there  exists  a  finite  M  >  N  and  the  sequence  of  vectors  [x  J  ^ 
has  N  linearly  independent  element  vectors  for  all  j  >  2M.  "...  If  we 
look  at  any  group  of  M  input  vectors  we  are  guaranteed  that 

||w  -Gii  <  jlw  -  G jj .  Since  M  is  finite,  this  occurs  an  infinite 

[I— j+M  —I!  |l— j  —II 

number  of  times  and  we  have 


1 im  w ,  =  G 

•j  _+  00  ^ 
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D.  CONVERGENCE  OF  THE  ADDITIVE-OUTFUT-NOISK  SYSTEM 

It  is  hard  to  imagine  a  real-world  process  that  does  not  have  some 
form  of  additive  output  noise.  In  addition,  one  inevitably  introduces 
noise  when  trying  to  make  measurements.  The  measurement  noise  can  be 
reduced  but  it  is  never  zero.  For  the  case  of  additive  output  noise 
n ,  we  have 

J  y  .  =  x  f  G  4  n  . 

.1  -  1 

and  thus  (2.2)  is 

e.  =  x  .  (G-w.)  +n 
J  — J  —  ~J  3 


Because  is  a  random  variable  it  is  necessary  to  look  at  the  behavior 

of  the  expected  value  of  the  weight  vector  w  as  i  _♦  »  .  To  start 

“•J 

with  let  us  consider  the  system  where  the  noise  process  is  zero  mean: 
Theorem  3.  Given  a  non- recursive  filter  G  of  length  N  with  additive 

■3 

output  noise  (lid,  zero  mean  and  finite  variance  o  )  and  an  N- 


dimensional  input  sequence  fx.]  ;  then  the  expected  value  of  the 

-jo 

weight  vector  w^  obtained  by  a-LMS  adaption  (0  <  a  <  ?)  con¬ 
verges  to  G  as  j  «  and  an  upper  bound  for  the  mean-square 
error  is 

_ _  2.T|,  .-,2 

o:  a  N  ix  i 

—max 


jlim^  [x‘  (w  .  -G)j2  < 


6(2  -  a6)  x  .  j 

— mmj 


Theorem  3  is  proved  in  Appendix  A. 


Thus  given  a  zero  mean  noise  source  which  is  uncorrelated  with  the 

.  OC  — 

input  sequence  { x ^ }  ,  the  expected  value  of  the  weight  vector  w. 

converges  to  G  as  j  -»  ™  .  In  addition  we  have  that  the  MSE  (mean 
square  error)  is  a  function  of  the  adaption  coefficient  i.  Thus  the 
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smaller  ct  is  the  smaller  the  MSE  will  be,  but  at  the  same  time  the 
slower  the  rate  of  convergence  will  be.  In  practice,  a  is  made  close 
to  one  during  periods  of  initialization  or  rapid  changes  in  the  physio¬ 
logical  process.  After  these  periods  a  is  then  reduced.  For  the 
experiments  described  in  Section  4,  ct  was  made  as  small  as  0.2.  Note 
also  that  the  MSE  is  a  function  of  the  maximum  and  minimum  input  vectors. 

The  fact  that  the  bound  is  divided  by  llx  .  I  is  not  unreasonable. 

||— min  !| 

Consider  the  example  where  we  have  a  zero  input  and  a  non-zero  output 
(due  to  the  noise  source).  In  this  situation  the  MSE  is  not  well  defined. 

Depending  on  the  system,  this  may  be  considered  an  undesirable 
feature.  In  addition  there  is  little  reason  to  believe  that  the  physio¬ 
logical  systems  one  may  want  to  model  have  zero  mean  additive  output 
noise.  To  handle  bias  in  the  observation  noise  it  is  necessary  to  intro¬ 
duce  a  non-linear  system. 

Definition.  The  augmented  weight  vector  w’  is  the  vector 


at  time  j,  where  the  augmented  input  vector  x'.  is 

1 

x  . 

-J 

The  new  weight  w  is  referred  to  as  the  bias  weight  because  under 
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suitable  conditions  its  expected  value  converges  to  the  mean 

value  (bias)  of  the  process.  Note  that  the  resulting  model  is  non¬ 
linear, 


N 

V 


y  ,  =  /  w  x  .  .  .  +  w 

3  ££  k  o 


Now  we  have  the  following  result: 


Theorem  4.  Given  a  non- recursive  filter  G  of  length  N  with  additive 

2 

output  noise  (iid,  mean  3  and  finite  variance  a  )  and  an  N- 


dimensional  input  sequence  {x^)^  then  the  expected  value  of  the 
augmented  weight  vector  obtained  by  a-LMS  adaption  (0<a<2) 

converges ,  namely  vi_  converges  to  G  and  converges  to  <3 

as  j  »  ,  and  an  upper  bound  for  the  mean  square  error  is 


lim 


x  <w  -  G) 


2  4  a  N  x 

<  II— max 


06(2  -o6) 


Theorem  4  is  proved  in  Appendix  A. 


The  augmented  a-LMS  Algorithm  thus  has  the  very  nice  property  that 

lim  w  =6  and  lim  w .  =  G.  Furthermore,  we  notice  that  the  norm 
o  — j  — 

J  “  J  “ 

squared  of  the  minimum  input  vector  is  one.  Thus,  for  all  finite  input 

2 

vectors  Xj ,  the  MSE  is  finite.  In  the  bound  for  the  MSE,  a  has 

been  replaced  by  its  upper  bound  of  4,  The  author  conjectures  that  the 

2 

MSE  bound  is  still  a  function  of  a  .  The  difficulty  of  working  with  a 
nonzero  mean  noise  arises  in  the  proof  because  the  a-LMS  Algorithm  is 
not  a  linear  process  ,  and  the  mean  3  is  not  known  a  priori. 
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E. 


CONVERGENCE  OF  THE  HIGH  AND  LOW  ORDER  MODELS 


From  both  a  practical  and  a  theoretical  point  ol  view  i t  is  impor¬ 
tant  to  know  how  the  OLMS  Algorithm  behaves  when  the  length  of  the 
model  is  not  equal  to  the  system  being  modeled.  Experimentally  it  is 
unlikely  that  the  length  of  the  impulse  response,  G,  will  be  known 
exactly  a  priori.  Theoretically  it  is  advantageous  to  make  N  as  small 
as  possible  because  this  reduces  the  MSE  bound.  The  usual  rule  of  thumb 
is  to  make  the  total  delay  time  of  the  tapped  delay  line  of  the  adaptive 
model  at  least  as  long  as  or  longer  than  G,  in  which  case  we  have: 


Corollary  1.  Given  a  non- recursive  filter  G  of  length  M  <  N  and  an 
N-dimensional  input  sequence  {x  )“  then  the  weight  vector 
obtained  by  a-LMS  adaption  (0  <  a  <  2)  converges  to  G  as 
j  _»  oo  where 


Proof  of  Corollary  1. 

Since  the  addition  of  feedforward  terms  of  zero  gain  does  not 
affect  the  output  of  G,  substitute  S  for  G  in  the  proof  of 
Theorem  2.  • 

This  completes  the  proof  of  Corollary  1. 

Thus,  when  the  model  is  of  higher  order  than  the  physiological 
system,  the  primary  effect  is  increased  MSE. 
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In  the  event  that  the  model  is  of  lower  order  than  the  physiological 


system,  the  problem  is  much  more  complicated.  Decompose  the  filter  G 


into  two  parts  G  and  G  where  G  is  the  first 
w 

and  G  is  the  remaining  elements  (see  Fig.  11). 

m  Tfti 

and  y^  =  xd^G  so  (2.2)  is  given  by 


N  elements  of 
In  this  way  y^ 


T  ~  t  »» 

e.  =  x  . (G  -  w  .)  +  xd .  G 

3  —  J  —  — j  — 3  — 


Now  unless  we  can  make  some  statements  anout  the  long  term  behavior  of 
the  input  sequence 


CD 


we  cannot  say  what  lim  w.  Is  equal  to.  To  see  this  we  make  the 
3  -»  00  ^  » s 

following  construction.  Suppose  we  interpret  y  .  as  being  an  additive 
output  noise,  then  y  must  have  a  stationary  mean  and  variance  if 
Theorem  4  is  to  apply.  While  this  may  be  the  case  for  some  systems, 
this  is  certainly  not  required  in  the  definition  of  an  N-dimensional 
input  sequence.  But  if  the  input  sequence  can  be  modeled  by  a  stochastic 
process  and  if  and  xd  .  are  uncorrelated,  then  Theorem  4  applies. 


This  being  the  case,  it  is  advisable  to  design  the  identification 
system  to  be  as  long  as  if  not  longer  than  the  maximum  possible  length 
of  the  impulse  response  of  the  physiological  system. 
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3.  FORWARD  TIME  CONTROLLER 


A.  THE  MINIMIZATION  PROBLEM 

The  problems  in  optimal  control  are  typically  associated  with 
dynamic  systems  evolving  in  time  where  the  systems  are  described  by  a 
set  of  differential  equations  [36-39],  Mathematically  formulated,  the 
determination  of  the  extrema  (optimal  points)  of  these  functionals  is 
a  problem  in  the  calculus  of  variations.  However,  in  the  case  of  the 
adaptive  model  controller  of  Fig.  7,  the  dynamic  system  is  modeled  by 


(3.0) 


This  mathematical  formulation  (parameterization)  of  the  dynamic 
system  in  conjunction  with  the  assumption  that  the  reference  input  r 
is  known  for  finite  future  time  leads  to  a  function  minimization  algo¬ 
rithm  which  does  not  use  the  calculus  of  variations.  The  development  of 
this  algorithm  is  the  thesis  of  this  section. 

Let  L  =  integer  constant 

N  =  the  number  of  weights 
j  =  the  current  time,  and 
T  =  LN  +  j 

The  reference  input  r,  and  a  set  of  constraints  a  ,  b,  ,  c  &  d 

i  iiii 

are  assumed  to  be  known  for  all  i  =  j+1,  j+2,  ...  T.  Thus  the  problem 
is  to  determine  x^  for  all  i  =  j+1,  j+2,  ...  T  such  that 

T 

2  (y.  -ri)2  (3.1) 

i=j+l 
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is  minimized  subject  to  the  constraints 


ai  < 


<  b. 


where 


c  .  <  x  .  <  d . 

1  -  i  -  1 


yi 


=  V 


k=l 


w,  x .  , 
k  i-k+1 


i  -  j+1 > j+2 ,  .  .  .  T 


The  blood  pressure  (output)constraints  a^  and  are  needed 

because  there  is  a  minimum  value  above  which  the  blood  pressure  must 
remain  at  all  times  (note  that  a^  <  0  has  no  meaning  in  this  problem! ) 
and  likewise,  there  is  a  maximum  value  of  blood  pressure  that  is  safe 
for  the  patient.  The  drug  rate  (input)  constraints  and  d  cor¬ 

respond  to  the  fact  that  there  is  a  minimum  and  maximum  drug  flow  rate. 
The  minimum  flow  rate  is  that  rate  which  will  keep  a  blood  clot  from 
forming  in  the  catheter.  If  a  blood  clot  forms,  the  catheter  must  be 
disconnected  from  the  drug  supply  and  flushed  to  clear  the  blood  clot. 
The  maximum  flow  rate  is  determined  by  the  maximum  amount  of  fluid  that 
can  be  added  to  the  blood  without  causing  undesired  side  effects. 

The  minimization  of  the  squared  control  error  (3.1)  is  used  for  two 
reasons.  First,  it  produces  very  satisfactory  control,  and  second,  it 
allows  a  considerable  reduction  in  the  mathematical  complexity  of  the 
forward  time  calculations  needed  to  achieve  closed  loop  control.  The 
forward  time  control  algorithm  described  here  is  thus  optimal  in  a 
minimum  squared  error  sense  with  respect  to  "controlling"  the  discrete 
non-recursive  filter  model.  The  optimality  of  this  approach  carries 
over  to  the  entire  control  system  when  the  model  has  converged. 
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Before  we  Begin  the  solution  of  the  minimization  problem,  it  is 


instructive  to  discuss  some  of  the  overall  system  requirements  and  model 
properties  which  have  been  used  in  choosing  the  most  desirable  algorithm. 
Notice  that  this  representation  of  the  physiological  system  (the  "plant11 
to  be  controlled)  gives  us  not  only  a  simple  formula  for  determining  the 
output  given  the  input,  but  it  also  gives  us  the  state*  of  the  model. 

This  is  very  important  in  achieving  control  because  we  are  dealing  with 
a  feedback  system  and  must  be  able  to  control  at  all  times.  That  is, 
we  do  not  want  to  wait  until  the  transients  from  the  last  correction 
have  died  out  before  we  make  the  next  correction. 

The  overall  system  requirements  arise  from  the  need  to  run  this 
model  reference  system  in  real-time  on  a  mini-computer.  These  require¬ 
ments  are: 

1.  The  algorithm  should  make  efficient  use  of  memory  since  the 
total  amount  of  memory  available  for  this  part  of  the  system 
is  small  (3-4K  words). 

2.  The  number  of  arithmetic  operations  should  be  minimized  be¬ 
cause  these  operations  are  not  very  fast  on  mini-computers  and 
there  is  a  limited  amount  of  time  available. 

3.  The  algorithm  should  calculate  the  solution  iteratively.  An 
algorithm  which  obtains  an  exact  solution  in  T  seconds  is 
useless  if  less  than  t  seconds  is  available,  while  an  algo¬ 
rithm  which  iteratively  refines  the  initial  guess  has  an  im¬ 
proved  "solution"  at  all  times. 


*The  model  is  a  non- recursive  digital  filter;  therefore  the  state  of  the 
model  is  precisely  x  . 
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4.  The  algorithm  should  converge  rapidly  and  in  a  finite  number 
of  steps  even  when  many  constraints  are  active. 

B.  THE  FORWARD  TIME  CALCULATION 

The  solution  of  (3.1)  is  divided  into  four  steps  (computation 

phases).  In  this  way  the  "solution"  of  (3.1)  is  refined  in  proportion 

to  the  amount  of  processing  time  available  at  time  j.  Physiological 

systems  quite  typically  have  some  amount  of  transport  delay.  This  means 

that  several  of  the  first  weights  along  the  non- recursive  digital  filter 

model  would  be  equal  to  zero.  In  practice  these  weights  are  not  exactly 

equal  to  zero,  but  instead  vary  by  some  small  amount  |  about  zero  (see 

Table  1).  This  is  handled  experimentally  by  labeling  all  |  <  £  as 

being  zero.  The  first  nonzero  weight  w^  is  then  referred  to  as  the 

first  weight  of  w  where  w.  now  has  N-J  elements.  At  the  same 
— J  ~3 

time,  the  model  is  iterated  forward  .'J  cycles.  In  this  way  the  model 
is  transformed  into  an  equivalent  system  with  no  delay, and  the  new 
is  not  equal  to  zero  (this  is  needed  because  we  will  have  to  divide  by 
to  get  the  solution  to  (3.1)). 

The  first  numerical  task  is  to  get  a  feasible*  solution  as  fast  as 
possible  (recall  that  the  amount  of  time  available  for  the  solution  of 
(3.1)  is  not  known);  This  is  done  by  recognizing  that  if  the  plant  is 
stable  then  it  will  have  a  steady  state  output  which  is  equal  to  r^. 

We  can  then  use  (3.0)  to  find  a  feasible  input.  Thus,  we  have 
Phase  1 : 

Find  the  steady  state  solution  x  and  then  let  x.  =  x  for  all 

ss  i  ss 

i  =  j+l,j+2,  ...  T.  Thus  we  have 

-fc* 

A  solution  is  feasible  if  the  constants  are  satisfied.  Thus  x^  is 

feasible  iff  c.  <  x.  <d,  and  a.  <y.  <b.. 

j -  i-i  i-i-i 
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LN 


x 

ss 


T 

y 

i=j  +  l 


ri 


N 

y 

k=l 


This  solution  is  interior  to  the  boundary  of  the  feasible  domain, 
but  it  is  not  likely  to  be  the  optimal  solution. 


Phas  e  2 : 

Check  to  see  if  the  problem  can  be  solved  with  zero  error.  Set 
2 

(y^  -  =  0  for  all  i  =  j+l,j+2,  ...  T  and  solve  for  x^: 


x. ,  c.  <  x.  <  d. 
l  l  —  l  —  i 


where 


undefined  otherwise 
N 

VJ  V 

-k+1 


-  /  W,  X  . 

i  .■H'  k  l- 
k=2 


w. 


(3.2) 


(  3.3  ) 


If  is  undefined  for  any  i  then  Phase  2  fails  and  we  proceed  with 

Phase  3;  otherwise  we  have  the  optimal  solution  for  (3.1)  and  are  done. 

Phase  3  is  a  heuristic  procedure  motivated  by  the  solutions  obtained 
using  Rosen's  gradient  projection  method  (Phase  4).  It  was  found  that  the 
solutions  for  the  optimal  drug  rates  obtained  in  Phase  4  could  be  fitted 
with  exponentially  decaying  envelopes.  Intuitively  this  is  reasonable 
because  we  would  expect  that  after  a  change  in  the  reference  input,  the 
drug  rate  would  change  greatly  at  first  and  then  settle  to  a  steady  state 
value.  This  heuristic  procedure  we  call  the  "clipped"  solution  because 
the  drug  rates  are  not  allowed  to  be  outside  of  this  envelope.  Thus, 

the  heuristic  procedure  is  that  if  the  drug  rate  3,  calculated  in  (3  31 
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is  outside  of  the  envelope,  it  is  replaced  by  the  value  of  the  envelope 
(see  (3.4)),  (Note  that  the  results  of  the  Phase  2  calculation  are  direct¬ 
ly  usable  in  Phase  3.) 

Phase  3: 


Check  to  see  if  the  "clipped"  solution  will  satisfy  the  output 
(blood  pressure)  constraints.  Solve  for  x^ ,  i  =  j+l,j+2,  ...  T  until 

y,  <  a.  or  y .  >  b. 
i  l  it 


where 


and  where 


L'B. 


x.  >  UB. 
1  —  1 


x.  =  /  x.  ,  LB  <  x.  CUB. 
l  \  i  l  l  l 


LB 


t’ 


x  <  LB. 
l  —  i 


A 


X  . 
3. 


r . 

i 


N 

y  .*  k T 

a  r  ‘■k*1 


(3.4) 


UB.  =  minfd. ,UC. ] 

l  L  l  iJ 

LB .  =  max [ c . , LC . ] 

l  L  i  l J 

UC.  =  x  +  T  expF-SHi-j  )  ] 

i  ss  L 

LC  .  =  x  -  T  exp[-'4<(i- j  )  ] 
l  ss  L  J 


where  B  is  defined  by 


*  =  NL  ln(T/^> 


where 


T=  max  x  -  x.  allowed  at  time  i  ^  j 
ss  i 
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a 


t  =  maxjx  -  x  allowed  at  time  a  -T 
ss  1 


If  Phase  3  is  successful  use  its  results  as  the  initial  guess  in  Phase  1; 
else  use  the  results  from  Phase  1.  In  order  to  describe  the  method  used 
in  Phase  4  it  is  necessary  to  formulate  (3.1)  as  a  nonlinear  programming 
problem.  Recall  that 

=  x[w.  i  =  j+l,j+2,  . , .  T 

or  in  matrix  form 


H 

CONSTANT 


where  the  past  inputs  x .  „  „  ...  x.  are  known  and  the  future  inputs 

j-N+2  j 

x.+1  ...  x^,  are  to  be  optimized.  Notice  that  the  upper  left  hand  block 
of  this  matrix  can  be  replaced  by  a  constant  vector  h  say;  thus  we  have 


and  defining  the  vectors,  we 

x  . 

-l 

,  X  <  j  +  N  (3.6a) 

,  x  >  j  +  N 

and 

for  m  >  k 

for  k  >  N  and  m  <  *  -  N  (3.6b) 
elsewhere 

Note  that  W  is  in  lower  triangular  band  form  with  bandwidth  N,  and 

that  all  of  the  elements  on  any  given  diagonal  are  identical.  Therefore 

using  (3.6b),  only  N  computer  words  are  required  to  store  W  instead 
2  ,, 

of  (LN)  words.  Furthermore,  inverse  exists  because  wi  is  not 
equal  to  zero  and  is  easily  obtained  by  one  front-solve*  operation.  We 
can  now  construct  the  following  quadratic  minimization  problem.  Minimize 


10. 


0 

0 


K+l-m 

,(at  time  j) 


or  letting  i  =  j  +  1  for  convenience, 
wri te  (3.5)  as 


w  he  re 


h  = 

K 


N 

Y 

|k=V-j 
0 


y,  =  h  +  if; 

—  1  —1  «Bt  * 


W.  ,  t  « 

k*l  j-k+1 


where 


V 

i=j  +  l 


(yi 


r.  ) 

l 


2 


yi 


=  h,  + 
l 


1 


10 .  .  ,  x 

i-j,k  j+k 


^When  the  matrix  of  coefficients  for  a  system  of  equations  is  lower 
triangular,  the  solution  is  readily  found  by  solving  the  equations 
one  at  a  time  from  the  top  down;  hence  the  name  front-sol’ e. 
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and 


=  reference  input 


a  <  y  <  b 
i  -  yi  -  i 


c  .  <  x  .  <  d 
i  —  i  —  i 


Define  the  error  vector  as 


z.  =  y.  -  r. 

-l  -i  — i 


=  li)x,  +  h .  -  r, 

=r—  1  -1  -i 


(3.7) 


Thus,  the  parameters  to  be  optimized  are 


x ^  =  h;  1  (i\  -  h^ )  +  W_ 


and  the  function  we  want  to  minimize  is 


'<£,>  -5<£A> 


(3.8) 


where  the  constraints  are  now  given  by 


a.  -  r.  <  zJ  <b.  -  r. 
—l  —l - i  -  —l  —l 


and 


c ,  -  _UI-1  <r„  -  h  )  <  JC_1z  <  d.  -  JT^r.  -  h.  ) 
— i  —  -i  -i  -  =  -i i  —  —l  —l 


(3.9a) 


(3.9b) 


The  desired  form  for  the  constraint  equation  in  the  nonlinear  programming 
problem  is 

-k-i  -  £k  k  =  1,2,  . . .  4LN 


33 


where 


-k-k  =  1  k  =  1 ,2 ,  . . .  4LN 

Rearranging  the  constraint  equations  (3.9a)  and  (3. 9b)  into  this  form 
we  get (without  normalizing) 


I 

r zn  > 

a .  -  r . 

-i  ~ 

—i  —I 

-I 

L> 

-  b .  +  r . 

—l  —l 

-  -  -  - 

JiT1 

c±  -  Jj)_  1  (r  -  h  ) 

-  -  -  - 

-A'1 

~d  +Li)1(r  -  h  ) 

—  _ 

—i  —  — i  — 1 

The  motivation  for  this  formulation  is  now  apparent  because  we  have 
that  the  gradient  of  f  is  given  by 

g  =  grad(f(zi))  =  z 

and  the  Hessian  matrix  of  second  partial  derivatives  of  f  is  given  by 

G(z.)  -  I  .  (3.11) 

One  can  now  apply  Rosen's  gradient  projection  method.  In  general, 
this  method  gives  linear  convergence  since  the  direction  of  steepest 
descent  is  taken.  However,  in  our  case,  where  the  Hessian  is  I,  the 
convergence  is  quadratic,  and  hence  the  minimum  of  f(z^)  in  any  given 
subspace  can  be  found  in  one  step.  A  general  description  of  Rosen's  algo¬ 
rithm  for  the  special  oase  where  G  is  the  identity  matrix, is  given  in 
Appendix  B.  A  programmed  version  of  Phase  4  is  given  in  Appendix  C, 
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Phase  4 : 


1.  Using  (3.7)  we  set 


where  z  denotes  the  error  vector  at  the  k  iteration  in  solv 
ing  for  the  optimal  z^,  and  is  given  by  (3.6a).  It  is 

known  from  previous  calculations  in  Phases  1,  2,  &3,  that  z° 
is  feasible  (all  constraints  are  satisfied)  thus,  the  number  of 
active  constraints  q  is  initially  equal  to  zero.  Also  compute 
U)  1  and  f  which  is  equal  to  the  right  hand  side  of  (3.10). 

As  a  matter  of  notation,  we  will  let 

=q  =  ££l’£2 . V 

T 

and  P  =  Q N  where  Q  Q  =  1  and  P  is  upper  triangular. 

=q  =  =q  =  =  — 

2.  Determine 


„  T  k 

i  .  -  n  .  z 

,  j  n  — j  n- 

.  .(active  I 

A  =  min  = 

\jn  -  T 

^  *  (constraints  1 

A.  >0 
jn- 

n  .  s 

-jn- 

3.  If  A*  >  1  then  go  to  step  6,  else  add  the  corresponding 


column  n 


to  N  and  form 


where 


and 


T  T 

P  7  =  N  n 
=q  -  -jn 


T  V 

P  =  (1  -  7  7)  2 


Set 


'nd 


q  =  q  +  1 

k+1  k  _* 
s  =  z  +  A" 


4.  Compute  the  vector  Q  where 


and 


Let 


T 

p  e  =  i 
=q  •£•  — q 


Iqfl  "  £ 


Cl  *  min 


.g)active  | 
J  (constraints) 


5.  If  (j  >  0,  then  let  s  =  -z  +  S  and  go  to  step  < 
delete  the  corresponding  column  from  ,  update  , 
q  =  q  -  1  and  go  to  step  4. 

6.  The  optimal  input  is  now  given  by 


c.  -  —  +  b  +  r  -  h  i 

i  wx  L  1  1  1 


If  x^  is  desired  it  is  given  by 


,,-l  ,  k 

xi  =  _Lj2_  (z  +  s  +  r  ~  h) 


This  completes  Phase  4, 


,  else 
set 
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By  making  use  of  the  speci  1  forms  that  some  of  the  vectors  and 
matrices  have,  we  have  th.  the  storage  requirements  are  as  follows 


QUANTITY 


MEMORY  LOCATIONS 
N 

2LN 


N  LN 

=q 

LN  (  LN  +  1 ) 

£  2 
)  NONE 


n  NONE 

-Jn 
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4.  EXPERIMENTAL  RESULTS 


A.  OPEN  LOOP  RESPONSE 

The  first  step  in  implementing  an  adaptive  model  controller  is  to 
obtain  open  loop  response  data.  This  is  necessary  to  find  out  if  the 
physiological  system  satisfies 

y(t)  =  F(x,y,t)  +  n(t) 

where  F(-)  defines  a  linear  deterministic  system  and  n(-)  represents 
a  noise  source.  If  F(-)  exists,  then  the  adaptive  model  controller  can 
be  used.  We  determine  N  and  A  (see  Section  2.B)  by  measuring  the 
open  loop  step  response  and  then  solving  for  the  impulse  response.  In 
our  experiments  on  dogs  this  is  accomplished  by  the  use  of  a  vasodepres¬ 
sor  agent,  Arfonad  (trimethaphan  camsylate)  which  produces  a  controlled 
state  of  hypotension  (shock).  Arfonad  is  primarily  a  ganglionic  blocking 
agent  which  has  the  effect  of  opening  the  dog's  blood  pressure  control 
loop.  When  the  dog  is  in  this  condition,  the  effect  of  a  vasopressor 
agent  such  as  Levophed  (1-norepinephrine)  is  easily  measured.  Figure  12 
shows  the  effect  of  a  step  change  in  Levophed  on  the  average  blood 
pressure.  The  two  step  responses  in  Fig.  12  were  taken  sequentially  on 
the  same  dog  with  the  time  between  steps  being  about  twenty  minutes. 

Using  this  data  it  was  decided  that  we  should  set  N  .-s  20  and  A  =  5 
seconds.  These  values  for  N  and  A  have  been  satisfactory  for  all 
the  dogs  that  we  have  run  on  this  system. 

Since  we  know  that  w  will  converge  to  the  impulse  response  G, 
we  would  like  to  know  what  G  would  give  these  ste?p  responses.  Let 
I^(t)  be  the  input  and  O^Ct)  be  the  output;  then  if  the  input  step 
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The  result  of  applying  (4.1)  to  the  step- response  data  shown  in  Fig.  12 
is  plotted  in  Fig.  13.  A  dominant  feature  of  this  data  is  the  transport 
delay  which  is  about  fifteen  seconds.  This  G  or  the  weight  vector 
from  the  previous  experiment  is  used  as  the  initial  weight  vector  W 
when  starting  a  new  experiment.  In  practice,  the  identification  process 
is  run  open  loop  until  a  good  model  has  been  formed.  When  the  MSE  is 
small,  the  loop  can  then  be  closed  for  automatic  control. 


B.  CLOSED- LOOP  RESPONSE 

The  instantaneous  blood  pressure  was  sampled  every  100  milliseconds 
and  then  averaged  over  a  30-second  period.  The  adaption  and  control 
cycles  occurred  every  5  seconds.  The  first  30  to  45  minutes  of  the 
experiment  were  typically  used  to  calibrate  the  data  link  and  to  perform 
open  loop  control  to  verify  that  everything  was  functioning  properly. 
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TIME  (SEC) 


Fig.  13.  The  Impulse  Responses  Derived 
from  the  Step  Response  Data. 

The  evolution  of  the  weight  vector  with  time  is  shown  in  Fig.  14 
(the  normalized  weight  values  are  tabulated  in  Table  1).  The  successive 
weight  vectors  are  plotted  off  to  the  right  to  aid  in  visualizing  the 
sequence  of  weight  vectors  as  a  continuous  three-dimensional  surface. 
Note  the  similarity  between  this  data  and  the  data  for  G  plotted  in 
Fig.  13.  The  times  at  which  these  weight  vectors  were  recorded  are 

indicated  in  Fig.  15a-g  as  "+"  marks  on  the  line  labeled 

_ 

(control  error)  . 


The  input  (drug  rate)  and  the  output  (average  blood  pressure)  are 
shown  in  Fig.  15a-g.  In  addition,  the  model  output  at  time  j  +  1  (one 
step  into  the  future)  is  plotted  as  a  cross-hatched  line.  The  line 


41 
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2 

labeled  (control  error) 
where  the  reference  input 
dashed  line. 


is  equal  to  the  average  value  of  (y  -  r  ) 


r  (pressure  set  point)  is  indicated  by  a 


In  Fig.  15a  &  b  one  can  see  not  only  that  the  system  regulates  the 
dog's  blood  pressure,  but  also  the  effects  of  learning  due  to  the  step 
changes  in  the  reference  input.  Now,  in  order  to  dramatize  the  inherent 
stability  of  the  system  due  to  the  continuous  system  identification  that 
is  occurring,  several  transients  have  been  introduced.  In  Fig.  15c  &  d 
the  effect  of  injecting  sodium  pentobarbital  (a  general  anesthetic)  is 
shown.  This  can  be  considered  a  typical  situation  because  patients  often 
require  drugs  in  addition  to  the  ones  used  by  the  control  system.  In  Fig. 
15e  &  f  the  dog's  blood  pressure  regulation  is  again  modified  by  the  ad¬ 
ministration  of  the  gasses  Halothane  (bromochlorotrifluoroethane)  and 
then  Amyl  Nitrite  (isoamyl  nitrite).  Halothane  is  an  inhalation  anes¬ 
thetic  and  Amyl  Nitrite  is  a  coronary  vasodilator.  The  reason  for  doing 
this  is  to  demonstrate  the  speed  at  which  the  system  can  compensate  for 
drastic  short-term  changes.  A  very  similar  situation  occurs  in  the  nor¬ 
mal  course  of  events  when  the  present  bottle  of  Levophed  is  running  low 
and  must  be  replaced  by  another.  The  time  during  which  no  drug  is  avail¬ 
able  can  be  minimized,  but  there  is  usually  a  detectable  difference  in 
concentration  of  Levophed  between  the  two  bottles.  Fig.  15g  shows  how  well 
this  system  regulates  when  the  physiological  system  is  not  undergoing 
drastic  changes. 

While  the  analytical  problems  introduced  by  a  time-varying  and/or 
nonlinear  system  are  still  unanswered,  the  empirical  results  shown  in 
Fig.  15a-g  are  a  strong  motivation  to  continue  with  this  approach  for 
regulating  physiological  systems. 
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Fig.  15a.  The  X-Y  Recorder  Output  Generated  During  the  Experiment. 
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Fig.  15b.  The  X-Y  Recorder  Output  Generated  During  the  Experiment. 
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The  X-Y  Recorder  Output  Generated  During  the  Experiment. 
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Kip.  15c.  The  X-Y  Recorder  Output.  Generated  IXirinp  the  Experiment. 
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Fig.  15f.  The  X-Y  Recorder  Output  Generated  During  the  Experiment, 
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The  X-Y  Recorder  Output  Generated  During  the  Experiment. 


CONCLUSIONS 


5. 


A.  SUMMARY  OF  RESULTS 

An  original  system  (the  adaptive  model  controller)  for  realizing 
an  automatic  therapeutic  control  system  has  been  presented.  It  was 
found  that  this  controller  can  be  realized  with  a  mini-computer,  and 
that  the  resulting  control  is  highly  satisfactory  from  a  medical  point 
of  view  [49], 

The  convergence  of  the  a-LMS  Algorithm  has  been  analyzed  for  deter¬ 
ministic  inputs.  Necessary  and  sufficient  conditions  for  unbiased 
convergence  were  given.  These  conditions  were  found  to  be  easy  to 
implement  and  to  cause  little  or  no  deterioration  in  the  control  problem. 
In  addition,  the  behavior  of  the  adaptive  model  where  the  system  to  be 
modeled  is  of  higher  or  lower  order,  was  presented. 

The  forward-time  controller  which  gives  minimum-squared-error  with 
respect  to  the  model  was  described.  The  minimization  was  done  by  a  non¬ 
linear  programming  technique  based  on  Rosen's  gradient  projection  method. 
The  solution  to  this  part  of  the  forward-time  calculation  was  optimized 
so  that  the  amount  of  memory  and  computer  time  required  was  at  a  minimum. 

An  actual  experiment  using  this  system  was  included  and  described. 
The  results  of  this  and  similar  experiments  have  been  very  successful. 

B.  RECOMMENDATIONS  FOR  FURTHER  WORK 

The  use  of  an  adaptive  multi-model  controller  (see  Fig.  16)  for 
blood  pressure  regulation  would  be  very  useful.  The  additional  adaptive 
processes  present  no  further  complication  in  the  identification  part  of 
the  system.  However,  the  optimal  controller  would  have  to  be  modified 
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Fig.  16.  Adaptive  Multi-Model  Controller. 


in  order  to  make  use  of  more  than  one  adaptive  model.  One  approach  would 
be  to  use  the  moctel  that  gives  the  "bestM  inpul-output  relation  for  the 
control,  and  to  use  the  remaining  models  as  constraints  only.  Thus, 
the  physician  could  not  only  specify  what  tin-  blood  pressure  should  be, 
but  he  could  also  specify  what  range  of  values  other  body  processes 
could  take  on. 

An  alternate  approach  would  be  to  extend  the  work  of  Ins lalle  ]  to 
decide  when  a  model  should  be  added  or  uroppeil  from  the  control  luop. 

A  controller  that  looks  very  promising  for  handling  nonlinear 
systems  is  the  adaptive  differential-model  controller  (see  Fig.  17). 

This  controller  was  proposed  by  Strom  ['18 J  who  conjectured  that  it 
should  provido  superior  performance  because  it  models  the  ’’local" 
behavior  of  the  system.  The  weakness  of  this  approach  is  thut  wo  must 
look  at  ’’derivatives”  to  form  the  model.  In  practice,  a  combination  of 
the  adaptive  model  controller  and  the  aduptive  differential-model  con¬ 
troller  may  prove  the  most  useful.  The  first  would  provide  a  "global” 
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Fig.  17.  Adaptive  Differential-Model  Controller. 

model,  while  the  second  would  give  a  " local "  model.  The  model  that  is 
used  to  find  the  optimum  input  would  then  be  a  function  of  the  magnitude 
of  input. 

The  inverse-model  controller  (see  Fig.  5)  and/or  the  feedback  model 
controller  (see  Fig.  6)  might  also  be  useful  for  controlling  physiologi¬ 
cal  systems.  Comparatively  little  is  known  about  these  controllers,  and 
hence  more  work  in  this  area  is  certainly  needed. 
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APPENDIX  A 


Norms  of  Vectors  and  Matrices 

The  notion  of  the  distance  between  a  vector  and  the  origin  is  called 
a  vector  norm  and  satisfies  the  following  properties  [2d, 45]: 

(1)  ||x|[  >  0,  ||x|[  =  0  iff  x  =  0 

(2)  ||cx||  =  jc  [  ||x||,  where  c  is  any  real  number 

(3)  ||x  +  y ||  <  ||x||  +  ||y ||  . 

We  will  use  the  Euclidean  norm  for  vectors 


IN  1 

I  4 

i=l  ii 


i  1/2 


Matrix  norms  have  an  additional  property  which  is 
(4)  ||AB||  <  ||A||  ||B||  . 

A  matrix  norm  is  said  to  be  compatible  with  a  vector  norm  if  for  any 
vector  x 


Ax  <  A  x 


(A. 00) 


The  Euclidean  norm  for  matrices  is 


A  = 


m  n 

I  I 

i=l  >1 


a/2 


which  is  compatible  with  the  vector  Euclidean  norm.  We  will  make  use  of 
the  fact  that  if  A  is  a  symmetric  matrix,  then 


||a||  =  max  |A,  (A) 


(A.  0) 
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Proofs  of  Theurfsi?  l,  2,  3  i.  4 , 


Proof  of  Theorem  1 , 


Given 


where 


-j  +  1 


w 


+ 


X  .€  . 
-.1  J 


e  .  =  xT(g  -  w .) 
J  -j  -  -J 


subtract  G  from  both  sides 


(w  -  G)  =  (w  .  -  G) 

-j+t  -  -j  - 


premultiplying  by  the  transpose  of  the  above  we  have 


(w  .  -  G)T(w  .  -  G)  + 
-J  -  ~ 


a 


(w  .  -  G)TX  .£  , 
~3  ~~3  3 


+  a 


2 

+  a 


T 

e  .x  .x  .e 

J-J-3 


3 


Thus 


therefo re 


a(2  -  a) 

>  0 


2 
e . 
j 


>  o 


liw  -  G  I  <  w  .  -  G  I 

!:~j+i  - 1  --  -j  - 1 


where  equality  holds  iff  =  0.1 

This  completes  the  proof  of  Theorem  1. 
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Proof  of  Theorem  2. 


By  Theorem  1  we  know  that  for  any  S-dlmensional  input  sequence 


tw 

fx,}  the  distance  between  w  and  G  is  non-increasing.  Geometri- 
-j  o  -j 


colly  this  means  that  we  can  construct  an  N-dimensional  sphere  in  E 


N 


centered  at  the  point  G  such  that  the  point  will  lie  inside  this 

sphere  for  all  i  >  j .  Thus  one  way  to  show  that  w  converges  to  G 

is  to  show  that  there  is  convergence  along  all  N  of  any  set  of  orthog 
N 

onal  axes  of  E  .  To  accomplish  this  we  introduce  the  operator  ( ■ ) 
that  operates  on  the  input  sequence  { x^ }  in  groups  of  N  to  produce 
sequences  of  orthogonal  inputs  {tj.J,  Thus,  for  j  =  2N,3N,  ...  let 


Jlj-N  =  -j-N 


7}.  =  (x.  e.)  e.  i  =  j-N+l,j-N+2f 

— »i  —  i  —i  —1 

where 

.  Jf  is, 

_1  i=J-N  4  \  '  •Di 


Si 


j-1  T 

l—i  "  .  vMTTil i 

1=J-N  -1! 


■  j-1 


L.(-)  performs  the  Gram-Schmidt  orthogonalization  on  our  linearly  inde 
J 

4  %t 

pendent  input  sequence  Thus  we  can  write  x_.  =  +  Xj  where 

tj.  is  the  new  information  and  is  orthogonal  to  x..  Thus 
J  J 


(w  .  -  G)  = 

-j+1  - 


i  -  a 


3iili 


<-  ~T 
x  .x  . 
J-'J  ,  „  J— 3 

i,2  a  ,2 

ft  N 


(w.  -G) 

-J 
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Premultiplying  by  the  transpose  of  the  above  we  have 


non-negative  definite 
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Now  looking  at 


2  a  and 


a 


we  see  that 


2  a  >  a 


for  o  <  a  <  2 


Thus  the  second  term  in  (A.i)  is  always  >  the  last  term.  Therefore 


r  T  " 

I  -a  — 

(w.  -  G) 

M  J 

"“J 

where  equality  holds  when  x^.  =  0.  The  convergence  due  to  the  {jjjJ's 
will  thus  form  the  desired  upper  bound.  To  find  this  upper  bound  ve 
need  to  look  at  the  convergence  due  to  after  N  adaptions. 

Thus  keeping  tabs  on  only  the  convergence  due  to  the  new  information  at 
each  adaption,  we  have 


I  -a 


(Wj-G) 


(A. 2) 


(“j+2  ”  — ^  = 


i  -  a 


-  2j+i 


where 


T  * 

i->*^ 

II  v  II2 

r  t  i 

i  -  a 

llv  (r 

1 

(w  .  -  G) 

-J  “ 

lhJ+lll  . 

IrJll  . 

-J+l 


0. 

-l 


4 


[terms  involving  x^ 


as  in  (A.I)] 
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(W  .  „  -  c)  = 

-J  +2  - 


I  -aj 


2  ,  ,9 

X.  A  lx  l 

-J 


!j«r 


(w.  -  G)  -  O.  , 
-J  “  0+1 


(w  -  G)  = 
-J+N  - 


I  -  a 


j+N-l  T 

V  Hi^i 


i=j  x 


(w.  -G)  -  0.  „  , 
-J  “  -J+N-l 


Thus  for  all  j=N,2N,3N,  ...  we  have 


(w  .  -  G)  = 
— J  — 


I  -  a 


o-i  t 

V  ii£ L 


i=j-N  ;x. 


%-*-&  -  2j-! 


Premultiplying  by  the  transpose  of  the  above  we  have 


w  . 

-0 


j-1  _T  ' 

V  nZi 
i  -  a  / - 

-  2j-l 

i=T-N  i[x±H2_ 

(A. 3) 


Now  given  that  0  <  a  <  2,  we  have  by  iteration  of  the  results  of  (A.  1) 


j-1  T 

\  HiHi 

(w ,  -G) 

I  -  a  /  - - 

,  2 

-0-N  - 

and  by  using  (A. 00),  we  have 


j-1  T 

V  ^i 

-a  Z  - 2 

1=J'N  INI 


W.  -  G 
-j-N  - 
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I  -  a 


j-l  T 

^  ii £l 


vsf  ^  iiiifrj-»-2ir  3  -2n’3k’ 


Now  to  find  out  what  ||Aj||  is  equal  to,  we  solve  for  the  eigenvalues  of 
A.,  Assume  that  c.  is  an  eigenvector.  Then 


Thus  we  find  that 


f  j-1  T\ 

I  -  a  2.  - 5  £,  =  A  c 

i-J-K  11*4*1  / 


Hj 


and  since 


1  «  W 

=  1  “  Of  TT-T 


S<M<, 

x  . 

-J 


we  have  that 


0  <  a  <  2 


-1  <  Aj  <  1 


H  -  maX  -X 

bj  "  j-N<i<j-l  Aj 


-1  <  b  <  1  for  all  j  =  2N,3N,  . . . 
J 
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Now  wo  can  write 


Let 

then 

where  5  >  0 
since 

we  have  that 

Therefore 


il*N-2li2  <  bo  ll^'2li2 
'l-2N~-l|2  <  bN  5  bNb!  11*0  -2il2 

j 

!K-£l!2<  n  bi  ll^-Gj|2  j  =  2N,3N,  ...  (A. 5) 

J  i^o 

2  2 
b  =lim  sup  b. 
max  r  l 

b2  2 

max  <  (1  -  a5) 

by  definition  of  an  N-dimensional  input  sequence;  and 


lim 
n  -»co 


<1  -  as>2n 


o 


iim 
j  -»«> 


.1 

n 

=r 


.2 

b .  w 
l  — o 


=  0 


lim 

w 

J-»°°  -j 


=  G 


.  I 


This  completes  the  proof  of  Theorem  2. 
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Proof  of  Theorem  3. 


Given 


we  have  that 


y  .  =  x  .G  +  n  . 
J  -J-  J 


w  =  w,  + 
j+1  “J 


a 


-j 


/  T„  T  \ 

x  x.G  +  n-xw 
j  “J-J/ 


Subtracting  G  from  both  sides  we  have 


(x  x  .  \  x  .n  . 

I  -  a  — ■ j  (w  -  G)  -  a  % 

l*i  7  Nl 


(A. 6) 


Taking  the  expected  value  with  respect  to  the  noise  source  we  have 

(recall  that  n.  is  zero  mean) 

J 


-5 


Now,  substituting  w.  ,  for  w.  ,  and  w.  for  w.  in  (A. 2)  we  have 
-J+1  -J+l  -J  -J 


lim  —  a 

w.  =  G  .  ■ 

J  _»  0°  -J  - 


This  completes  the  first  part  of  the  proof. 


G3 


Let 


d.  =  I  -  a 


T 

X  X  . 

--J-3 


x 


Then  (A. 6)  is 


and 


Vj 


(w .  -  G)  =  D. <w .  -  G)  -  7. 

-j+1  -  -3  “3  ”  -3 


(A.  7) 


and  we  have  for  an  N  step 


(— j  +2  ~  — )  =  £j+l(^+l'S>  “  2j+l 


=  D.  .  (D  (w  -G)-  7.)  ~  7., , 
=3+1  =j  -j  -  -3  -3+1 


Sj*j5)<=j  -5’  - 


-2Jrf 

■  ?:j*2iViV-r£>  -  Sn-Jitii]  -  Jj+A-1-1  ■  2J*2 


j+N-i  v-1  N-i 

^t.s'5’  *  jl  2k '-j  "2>  -  j4  kj]_1 


Thus  for  all  j  =  N,2N,3N,  ,  we  have 


(w  .  -  G) 

-3  - 


j -1 


N-l 


k=j-N 


[I  ..  2^-N-V  ~  1  ,  J_.  =j-N+k— j-l-i  -  2j-l 


(A. 8) 
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Premultiplying  by  the  transpose  of  the  above 


2  jrr1  f2 

=  n 

J  k=j-N  " 


+  |j-J-l||2  +  1=^-!— J-2|j2  +  + 


X 


+ cross  terms  involving  7^ ' s 


Taking  the  expected  value  with  respect  to  the  noise  source,  the  cross 
terms  involving  y  drop  out  because  the  noise  is  zero  mean,  thus 


w  .  -  g|  = 
t-j  -I 


jn 

n 

k=j-N  J 


+  l|tj-if  *  Ift-ilj-sf  +  * 


kjL 


2  i  2  2 

+  M  +fe-i  M  ♦ 


k  ^2  =k| 


(A. 9) 


Now  making  use  of  (A. 0),  we  have 


£j-i||2’  ||=j-i=j-2||  •  ’ 


JnL 

k=j-N+2 

£k 

=  A  =1  (A. 10) 

max 
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Therefore 


k=j-N 


Noticing  that 


OCX-  ,n  . 
-3  3 


T 

2  2  -j-j 
a  n  .  — i-i. 
J  llx  H4 
i~j 


2  2 
a  a 


we  can  write 


n 

k=j-N 


2  2 
a  o 


N 

/ 

‘  -  - 

k=l 


Now  we  have 


/ff  &<2! 

k=j~N 


•j-N 


„  2  2 
Na  o 


^min 


Now  referring  to  (A. 3)  we  see  that 


j-1 

k=‘j-N  k  -1  N 


can  be  written  as 


>-c  f 

i=j-N 


T 

^i 


(Wj.N-G) 


Vl 


Thus,  given  that  0  <  a  <  2,  we  have 


w 

~3 


< 


+ 


2  2 
Na  0 
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where  A^  is  as  defined  in  (A.4).  Now  to  get  an  equation  similar  to 
(A. 5)  we  do  the  following 


where 


2  2  2 
w-G  <b  w-G  +w 

— N  —  —o-o  — 


2  2 

Na  a 


— ruin 


(A. 11) 


^2N"-|2  ^  bNK'£l|2  +  “ 


2  2  2  2 
<  b  b  w-G  +  (b_T  +  l)w 
—  No  — o  —  N 


— 3N  ~  -  -  b2N  —  2N  "  -  +  U 


<  b2NbNbo  +  (b2NbN  +  b2N  +  1)u 


|  w .  -  G  <  ff  b  w-G  + 

l-j  -  -  1  -°  - 


2 


i=i  k=j+i-i 


a  .<*• 


w  j  =  2N,3N, 


Now  substituting  b  for  b  we  get 
max  l 


|w .  —  G|| 2  <  b2J  ||w  -G||2  +  u  )  b2J 

• <-*  max 

i=o 


-||  -  max  ||  o  -|| 
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Now  since 


llm  b2J  =  0 
j  ->®  max 


anti 


lim 


T  „2‘  . 


1  -  b 


max 


we  have  substituting  for  w  and  taking  the  limit 


lim  i  2 

I  w  .  -  G  < 
J  -*  00  -J  "  ~ 


2  2„ 

a  o  n 


d-b2  )  |  x  2 
max  —min 


2 

ao  n 


5(2-a8) 


f^min 


and  using  the  Cauchy-Schwartz  inequality  we  have 


and  thus 


hllh-sl 

11"  |W-c>]2  < 

j  _>  oo|_  J  -J  J  5(2-  CX&)  jjx  . 


i  — min 


This  completes  the  proof  of  Theorem  3. 
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Proof  of  Theorem  4. 


Given 


ax\ 

w'  .  =  w'.  + 

~j  +  1  "J  ||  x' 

-j 


-J—  (xT  G  +  n.-x'Tw,'j 
,  ||  2  y-J  -  J  —3—3/ 


we  have 


T 

w',  ,  =  w'  +  - —  x'. 

-J+l  -J  xi  I!2  -3 
“j 


n .  -  w  . 
J  °J 


G  -  w . 
-J 


(A. 12) 


where  w  is  weight  zero  at  time  j;  or  in  partitioned  form 
oj 


W  ,  , 

OJ+1 

"■  w  * 

Oj 

J. 

a 

w  .  - 

L  ^ 

V/  , 

„  -J  . 

T 

i  2 

X  .  1  + 

.  “J| 

1  '  x  . 
.  “J 


T 


I  T 


Taking  the  expected  value  with  respect  to  the  noise  source  we  have 


w'.  ,  =  w'.  +  a 

-j+i  -j 


X  '  X  ' T 

-3-j 


L-j 

1  -3-3 

the 

noise 

r 

3 

-  w  . 

OJ 

G 

-  w . 

-J 

n  .  -  w 

J  OJ 


G  -  w  . 

L  ~  "J  J 


Now,  subtracting  G'  from  both  sides  where 


G'  = 


e 

G 


we  have 


(w T,  -  G'  )  = 
-J+l  - 


x'x; 


i  -  a 


(w  -  G’  ) 
J 


Now,  making  the  appropriate  substitutions  for  w,  G  and  x  in  (A. 2) 
we  have 

lin  IIS',  -  o' II2  =  o 


J  -»  °°I|-J 
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or  equivalently 


11m  . 

.  w  .  =  P 

j  -oj 


lim  —  _  . 

w  =  G  .  ■ 


This  completes  the  first  part  of  the  proof. 
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Letting 


x’x* 

d;  =  i  -  a 

INI 


n  ,  -  n  .  , 
j  J+l 


4  ’ 


0  J 


we  can  write  (A. 12)  as  follows 


=  w .  , 
-J+1 


aj  x’x’T 

-  a  — (w;  -  g'.)  +  r 
Q  I  ^  ^  ^  ^ 


(w'  ,  -  G'.  ,  )  =  D'  (w -  G’)  +  r 

-  j+i  -j+i  -i  -a  -  -J 


Now  substituting  into  (A. 7),  we  have  for  an  N  step 


j-1  ^  N-l 

II  +  2  f]  2i- 

J  J  k=j-N  J  3  i=l  k-N-i  J 


y  +  7  ’ 

N+k-j-l-i  -j-1 


for  all  j  =  N.2N.3N,  ... 
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I’reiiiu  1 1  ip  1  y i ng  by  the  transpose  of  the  above 


9 

+  (N  -N  cross  terms  involving  2^  ■  7^  <  i  ^  k)  ( A .  1  it  1 

Consider  individually  the  expected  value  with  respect  to  the  noise  source 
of  the  four  terms  on  the  right. 

1)  The  first  term  of  (A. 13)  is  (referring  to  (A. 3)) 


2)  Consider  the  second  term  of  (A. 13): 
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referring  to  (A. 9)  and  (A. 10),  we  have 


but  recall  that 


second  term 


5.2  iiz 


k=i 


-j-k| 


=  <ni-ni+i)2 


Therefore 


2  „  2  2 
ni  -  20  +  Vl 


=  2  a 


second  term  <  2No 


3)  Consider  the  third  term  of  (A. 13): 


iT/N-1 


Note  that  only  the  2j_ff  term  is  correlated  with  thus  all  the 

other  terms  are  equal  to  zero  because  J\  =  [0  0  ,  .  .  0]  .  Thus 


third  term  =  2(w 


”j.N-Sj-N>T  _k  Si  (I  S-lnUj-K 


or 


where 


2<w'_n-G’VTDD.7._n 


■j-1  *ITN-1 

°D-  =  n  o;  n  d-  « ,, 

— J  [kJ/.N  Lk=1  =J-N+k 
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In  expanded  form  we  hove 


third  term  =  2(n.  -n.  )(w’.  -  G '  „)  DD.Z 

j-N  J-N+1  -j-N  -j-N  =j- 


Z  =  [100  .  .  .0] 


The  w\  ..  vector  drops  out  because  it  is  uncorrelated  with  n  and 
-J-N  J-N 

n j-N+1  anc*  is  therefore  multiplied  by  zero.  Similarly,  only  the  first 

term  of  G'  .  „  in  non-zero.  Thus 
-  J-N 


third  term  =  -2 (n  .  ,  -  n  „  ,)(n.  „)  Z  DD.Z 

J-N  J-N+1  j-N  -  ==j- 


2  2  T 

-2  (n  .  -  &  )  Z  DD.Z 

j-N  -  — j- 


2  T 

-2  0  Z  DD.Z 

J 


Now,  making  use  of  (A. 00),  we  have 


2  T 

third  term  <  2o  |z  DD^.  Z  | 


<  2o  Z'j  ii  DD.Z 


<  2(?  iizll  jDDjl  j|Z|| 


(A.  14) 


liZl  =  1 


"  1 


therefore 


third  term  <  2o 
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4)  Consider  the  final  term  of  (A. 13): 


(N2 


-  N 


cross  lerms  involving  7',  7^,  i^k) 


Notice  that  these  cross  terms  can  be  written  in  the  form 


(n  -  n,  )(n  -  n,  )  Z  M2 
1  l+l  k  k+1  —  -» - 


where  M  is  a  matrix  corresponding  to  the  product  of  the  Dj's.  Only 
2(N-1)  of  these  terms  are  non-zero.  These  are  the  terms  where 
|i-k|  =1.  For  each  of  these  terms  we  have  (see  (A.  14)  and  following) 


(n .  ~  n  .,  ) (ii  -  n  ,)  Z  MZ  = 
1  1+1  k  k+1  —  =  — 


2  T 

-  0  Z  MZ 


<  c“  Z  Ml!  Z 


<  c 


Therefore 


We  have  now  proven  that 


fourth  term  <  2(N-l)o 


2  2  .22 
w'-G'.  <  A'.  w’  -G'.J  +  4o  N 

-J  J  ~  =J  -J-N  —  J-N  | 


for  all  j  =  N,2N,3N,  . . . 


Now  substituting  into  (A. 11),  we  have  the  desired  result 


lim 


T  "I2 

x;  (w'_  -G’)  < 

|_“J  -J  -3 


2  4o  N  x’ 

— max 


o:5(2  -  ctB) 


This  completes  the  proof  of  Theorem  4. 
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APPENDIX  B 


A  General  Description  of  the  Algorithm  Used  in  Phase  4. 

In  the  last  decade  many  methods  have  been  developed  for  minimizing 
the  quadratic  objective  function  (3.8)  subject  to  linea.  inequality  con¬ 
straints  (3.10)  [40-44].  One  of  the  most  promising  is  given  by  Goldfarb. 
His  algorithm  extends  Davidon’s  variable  metric  technique  for  constrained 
minimization  and  is  based  on  ideas  found  in  Rosen  [41  ].  When  the  Hessian 
is  I,  as  it  is  in  our  case,  Goldfarb' s  algorithm  is  equivalent  to 
Rosen's  Gradient  Projection  Algorithm. 

The  algorithm  presented  here  is  specifically  designed  for  the  case 
of  a  quadratic  objective  function  (3.8),  linear  inequality  constraints 
(3.10),  the  Hessian  matrix  of  second  partial  derivatives  being  equal  to 
the  identity  matrix  (3.11)  and  the  gradient  being  equal  to  z. 

The  algorithm  described  in  Phase  4,  which  incorporates  the  ideas  of 
Rosen  and  Gill  and  Murray  [47],  was  worked  out  by  Kaufman  [46]. 

We  want  to  minimize 

k  ifT  k 

f(£>  =  4(zk  zk) 

subject  to  the  constraints 

nTzk  >  l  m  =  1,2,  ...  4LN 

— m—  —  m 

where 

nTn  =  1  m  =  1,2,  .  .  .  4LN 
— m— m 

A  necessary  and  sufficient  condition  for  f(z*)  to  be  the  global  minimum 
is  that  there  exists  an  (3  such  that 

g(z»)  = 
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where 


N 

=q 


■2’ 


n  } 

-q 


The  columns  of  N  are  the  q  unit  normals 
=Q 

whose  intersection  is  an  affine  subspace  in  which 
tion  equation  is  thus 


where 


P  =  I  -  N  (NTN  )'  *N" 
-q  =q  -q-q  -q 


to  the  hyperplanes 
z*  lies.  The  itera- 


(B.  1) 

CB.2) 


is  Rosen's  projection  operator  that  projects  Em  into  the  constraint 
manifold.  Note  that  in  the  unconstrained  case  q  equals  zero  and  P  =  I 

-q 

k+1 

and  thus  (B.l)  reduces  to  the  one-step  Newton  method,  z  =0. 

Substituting  for  P  in  (B.l)  we  have 
s=q 


or  equivalently 


k+1 

z 


T  -1  T  k 
N  (N*N  )  N  Z 
=q  =q=q  <-q— 


k+1 

z 


n  a 

=q- 


(B.3) 


where 


a  =  (N>) 

-  =q»q 


(B.  4  ) 


The  elements  of  correspond  to  the  q  active  constraints  of  t. 

The  computation  of  (2  can  be  quickly  performed  as  two  back  solve  opera- 

T 

tions  if  the  following  construction  is  employed.  Assume  that  is 

T  T 

q  X  t,  then  N  can  be  written  as  N  =  QP  where  Q  Q  =  QQ  =1  and 

=q  =q  =q  ■»■**=  = 

is  a  q  X  q  upper  triangular  matrix  [45].  Thus,  (B.4)  can  be 


written  as 
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N  N  a  -  2 

=q=q—  —q 


T  T 

£qg  2lq«  =  £q 


p  p  a  -  £ 
=q=q-  — q 


Since  is  an  upper  triangular  matrix,  _  can  be  quickly  calcu¬ 

lated  by  using  two  back  solve  operations. 


PrS  -  £ 
=q-  -q 

and 


The  addition  of  constraints  to  N  is  done  as  follows; 

«q 


Let  N  , 
=q+l 


*  n  1 

L-"»  ,  -qj 


If  N  =  QP  ,  then 
-q  =  -q 


T 

N*  N  , 
»q+l=q+l 


- 

- 

.  T 

» 

T 

N  N 

N  n 

=q-q 

1 

=q-q 

T 

» 

T 

n  N 

n  n 

-q=q 

1 

_q_q 

l 

-* 

T 

1 

,T 

P  P 

N  n 

=q=q 

1 

=q-q 

T„ 

1 

T 

n  N 

n  n 

-q=q 

1 

-q-q 

- 

1 

_ 

If  we  put  Pq+j  in  uPPet  triangular  form,  then 


lq+1 


P  '  7 

-q  ,  - 


0  i  P 
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and 


T  T 

N  . N  =  P  „P  = 

=q+l=q+l  =q+l=q+l 

P  P 
=q=q 

!  PT  7  " 

1  =q  - 

T 

7  P 

-  «q 

1  T  2 

i  j  y  +  p 

Thus  to  calculate  P  ,  we  need  to  solve  for 
=q+l 


y  and  p  where 


and 


T  T 

P  y  =  N  n 
=qi.  -q-q 


=  (n^n  -  7 T7)^ 
-q-q - 


=  (1 


T  4 
7  7)2 


When  a  constraint  is  dropped  the  corresponding  column  must  be 

th 

deleted  from  P  .  If  the  k  column  is  deleted,  we  have  q  -  k  ele- 

=q 

ments  below  the  diagonal  which  must  be  set  to  zero.  This  is  easily 
done  using  Given's  rotation  matrices  [45]. 

The  author  suggests  that  interested  readers  see  [40 ]  Section  8 
and  [46]  for  a  discussion  of  rates  of  convergence  and  operation  counts. 


APPENDIX  C 


A  PROGRAMMED  VERSION  OF  PHASE  FOUR 


real  PROCEDURE  PHASE4(IM,LN,W,X0,R,A,B,C,D)  f  VALUE  N,LNI 
INTEGER  N,LN!  PEAL  ARRAy  w,X0,R,A,B,C,D; 

BEGIN 

REAL  T1 ,T2,T3,T4,LSTAR,ASTAR,NORM2S,NQRM2ZF 
TNTEGER  0/ Jl,  J2, J3, 1 1,12, 13,  INDEX,  START! 

8  PDIM0££LN(LN*n)/2J 

FOUATL  MAXLN943, LN401 1>0,  X V*?l ,  0 ' -2, PDIHP820, EPS ILONO0 . 25  * -5 1 
REAL  ARRAY  HR, «IN,Z,S, NORM, ALPHA, aETAlUHAXLN]  ,PtllPOlMJ  , 

L  !  1 !  L  N  4  J ; 

INTEGER  ARRAY  MQ  1 1 j HAXLN]  ) 

INTEGER  PROCEDURE  IN  1 ( I )  J  INTEGER  It  INIOC  C 1-1 3  MOD  LN)+i; 
INTEGER  PROCEDURE  IN2(I);  INTEGER  It 
CASE  (((I-I)  0  LN5+1  ) 

BEGIN 

IN20-2;  IM2C-1J  IN202I  IN201; 

END) 

BOOLEAN  B1,B2; 

BOOLEAN  ARRAY  ACTIVE (1 S  L  N  4 1  / 

LAPEL  STEP! ,STEP2,STFP3,STEP4,STEP5,STEP6I 

8 

&  v  *  REPLACEMENT  0  *  INTEGER  DIVISION 
8 

8  WIN*  WB  INVERSE 
8  WIN(I,J)s  W(I*l«i)  FOR  1>«J 
&  P(I#J>«  P(I*J(J-l)/2)  FOR  I<*J 

8  INI  AND  IN2  ARE  INDICATOR  FUNCTIONS  SUCH  THAT 
8  IN2 ( J  3  *  -2  IF  CONSTRAINT  IS  I 

&  IN2CJ3*  -1  IF  CONSTRAINT  IS  -I 

8  IN2CJ5*  2  IF  CONSTRAINT  IS  WIN 

8  IN2(J)*  1  IF  CONSTRAINT  IS  -WIN 

8  IN  1  ( I )  IS  AN  INDICATOR  ARRAY  WHICH  EQUALS  THE  ROW  NUMBER 
8  OF  THE  CORRESPONDING  IN2CI)  BLOCK, 

8  NO (13  IS  an  INDICATOR  ARRAY  WHICH  POINTS  TO  THE  ITH 
8  ACTIVE  CONSTRAINT 

8  NORM  (I)*  NORM  OF  THE  I-TH  ROW  OF  WIN 
8  ACTIVE (I)*  TRUE  IF  THE  I-TH  CONSTRAINT  IS  ACTIVE 

8  INDEX*  THE  COLUMN  WHICH  IS  TO  BE  ADDED  OR  DROPPED 

8 

STEPi:  OOP;  START*?! ! 

&  COMPUTE  HRSH-P 
FOR  JVl  TO  (NlVN-l)  DO 
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BEGIN 
TiO-R  tJJ J 

FOR  KON  STEP  -1  UNTIL  (Jl©J+i)  DO 
TIOTUW  t K J  *X0  [J1-K3  | 

hr  cjjoti; 

END  j; 

FOR  JON- TO  LN  DO  HR  C J3  t J1  F 

&  COMPUTE  Z0*  W8*X0*H-R 
FOR  JOl  TO  LN  DO 
BEGIN 

TlOHR  t J] i  Ji©J*i; 

FOR  KOCIF  J<«N  THEN  1  ELSE  J+l-NJ  TO  J  Du 
T19Tl+WtJi-KJ*X2£KJ; 

ZtJlOTl/  SlJJO-TlJ 
END  J/ 

&  COMPUTE  WIN  &  NORM 

NORM  til ©T30WIN tll©t,0/W  111  J  T2©T3*T3/ 

FOR  J©2  TO  LN  DO 
BEGIN 

T1©0,0>  JIOJ  +  H 

FOR  K©(  IF  J<«N  THEN  l  ELSE  J*1»N)  TO  (J20J-1)  DO 
T19TWWIN  IK)  *W  [Jl-Kl  I 

T40-T3*Tll  wINtJJ©T4|  T2©T2*T4*T4J  NORM  £ J1 ©SORT (T2) I 
END  J; 

&  CALCULATE  L 
FOR  J©1  TO  LN  DO 
BEGIN 

L  [  JJ  ©A { J] »R  [  J] |  L tJ*LN10R£JJ«BtJl )  T1O0.0J  J1©J*1> 

FOR  KOI  TO  J  DO 
TlOTUWINlJl-KMHRCKlI 
L  tJ+2*LN]OC  tJ]*Tl;  L  l  J+3*LN]  ©*(D  IJJ  *T1)  t 
FOR  K©0  TO  3  DO  ACTIVE £J+K*LNJ ©FALSE* 

END  J; 

STEP2I  LSTARO1.0I  Nl©4*LN* 

FOR  J©1  TO  Nl  DO  IF  NOT  ACTIVE  IJJ  THEN 
&  CALCULATE  LAMBDA [JJ  FOR  ALL  J  NOT  IN  THE  SET  OF  ACTIVE 
*,  CONSTRAINTS 

BEGIN 

JlOINl(J)f  B 10IN2 ( J 1  * 

1 1©IF  A8S(ai)sl  THEN  -l  ELSE  H 
IF  Bl  THEN  TlO(Il*LtJl-ZIJll  1/StJlJ 
ELSE 
BEGIN 

T3OT2O0.0*  J20J1+1* 

FOR  KOI  TO  Jl  DO 
BEGIN  T49WJN CJ2-K] ) 

T30T3*T4*Z £K] I  T20T2*T4«S tKJ I 
end  ki 

T 10( 1 1 *L  IJ1-T31/T2I 

endj 

IF  CT1>EPSIL0N)  AND  <T1<L<STAR3  THEN 
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&  RECORD  THE  VALUE  AND  INDEX  OF  THE  MINIMUM,  NON-* 

&  NEGATIVE  VALUE  OF  LAMBDA) 

BEGIN 

lstahoti;  INDEXOJ) 

END? 

END  J; 

STEP3J  IF  LSTAHal  THEN  GO  TO  STEPS) 

8,  ADD  THE  INDEX  JUST  FOUND  TO  THE  SET  OF  ACTIVE  CON- 

8,  STRAINTS 

QOG+1?  STARTOQ?  nq cg] sindexi 

IIWIF  ADS ClM2(INDEXi ) *1  THEN  -LN  ELSE  LN? 

ACTIVE  CINPEX  + 1  INACTIVE  t  INDEX]  <?TRUE) 

S,  UPDATE  2 

FOR  JOl  TO  LN  DO  Z C J3 <?Z f JJ +LSTAR*S UJ I 
IF  Qa  1  THEN  BEGIN  P11391)  GO  TO  STEP4)  END) 

&  UPDATE  P  BY  ADDING  A  COLUMN 

B1CIN2CINDEX] )  N 190-1 1  J3«?Q*N102J  1 39IN1 C INDEX ) )  T3C0.0I 
FDR  J«>1  TO  N 1  DO 
BEGIN 

1 19NQ  t J3  ;  I2<?IN1(I1))  JlOJ-1)  J20J* Ji 02)  B20IN2CI1)) 

IF  B2  THEN 
IF  B1  THEN  T1O0.0 

ELSE  Tl^IF  1 2> 1 3  THEN  0  ELSE  WIN tI3-I2  +  i] /NORM tI3J 

ELSE 

IF  B 1  THEN  T 40 IF  I3>12  THEN  0  ELSE  WIN  112-13*1] /NORM CI2] 
ELSE 
BEGIN 

T1T>0,0)  I10ARSCI2-I3) ) 

FOR  KOI  TO  (IF  1 2 < 1 3  THEN  12  ELSE  13)  DO 
TJOTI*WIN[K]  *WIN£K+IU) 

T 1CT  1/ (NORM £12] -NORM [13]  )) 

END) 

IF  ABS (B2*B 1 3  *2  THEN  Tl^-Tl) 

FOR  KOI  TO  J1  DO  T1^T1-P [K*J23 *P  IK*J33 ) 

T20P  CJ  +  J33  OT i/P  IJ*J2J ) 

T30T3*T2*T2 I 
END  J) 

P  IQ  +  J3]OSQRTC1,0-T3)? 

STEP4;  FOR  JOSTART  TO  0  DO 

&  SOLVE  FOR  THE  ELEMENTS  OF  BETA 
BEGIN 

1 1  ON Q [ J ] ?  J10J-i?  J29J*J102 )  T40LII1J) 

T 1 91 F  IN2CII)  THEN  T4  ELSE  T  4/NORM  (INl(Il)]) 

FOR  KOI  TO  J1  DO  TA$TI-P  IK  +  J2J *BETA  [K] ) 

BETA  [J]CTl/PtJ*J2]  ) 

END  J) 

ASTARO0.0) 

FOR  J9Q  STEP  -1  UNTIL  1  00 
&  CALCULATE  THE  ELEMENTS  OF  ALPHA 
BEGIN 

T10BETA [J] )  J1SJMI  J30J20J*JI02) 
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FOR  KOJ1  TO  0  DO 
BEGIN 

T10TI-.P  £J*J2J  *ALPHA  tK]  >  J29J2*K 

END> 

T10ALPHAUJ0T1/P[J3]  I 

&  FIND  THE  MXH'MUM  ELEMENT  OF  ALPHA  AND  REMEMBER  ITS 
&  INDEX 

IF  T1<ASTAR  THEN  BEGIN  ASTAROTiJ  INDEXOJJ  ENDl 
END  J; 

STEP5I  IF  ASTAR*0  THEN 
BEGIN 

NORM2SCNORM2Z90.0> 

FOR  JOi  TO  IN  DO 
&  CALCULATE  THE  ELEMENTS  OF  3 
BEGIN 

TIO-Z [ JI J  NORM2ZONORM2Z*Tl*Ttl 
FOR  KOI  TO  Q  DO 
BEGIN 

I30NQ  £ K 3 }  I20IN2CI3)!  Ii0INi(I3>» 

T20IF  12  THEN  IF  I !■ J  THEN  l  ELSE  0 

ELSE  IF  J<»It  THEN  WIN  IIUI-JJ  /NORM  II  t]  ELSE  01 

IF  ABS ( 12 ) ■ 1  THEN  T20-T2I 

TiOTl*T2*ALPHAJKI> 

END  k; 

S[J30T1>  NORM2SONORM2S*T1*TIJ 
END  J) 

IF  N0RM2S  <XV*NORM2Z  THEN  GO  TO  STEPS 
ELSE  GO  TO  STEP2I 
END  IF  I 

S  DROP  THE  APPROPRIATE  COLUMN  FROM  P 
QOQ-iJ  STARTOINDEX  }  I20NG (INDEX) I 
1 157IF  AB3(IN2(I2})«1  THEN  -LN  ELSE  LN» 

ACTIVE [12+11  INACTIVE [121 OFALSE I 
IF  INDEX«Q+1  THEN  GO  TO  STEP4 I 
FOR  JOINDEX  TO  Q  DO 
BEGIN 

JlOJ*lf  I10J2CJ1*J/2J  J30J2-JJ  NQtJJONQtJl)  I 
FOR  KOI  TO  INDEX  DO  P [K+J3J 9P  IK*J23 » 

&  USF  A  GIVENS  ROTATION  TO  ELIMINATE  THE  EXTRA  ELEMENT 
T  iOP  t J* J2J I  T20PIJ1+J2JI 
P  t J23  0T30SQRT (Ti*Tl*T2*T2)J 
T10TI/T3J  T20T2/T3/ 

FOR  KOJ1  TO  Q  DO 
BEGIN 

1201 1 )  !10II*K>  T40PlJ*m>  T39PCJi*Il}> 

P  CJ+I2JOTi*T4*T2*T3>  P  [Jl+Ill OT2*T4-Tl*T3l 
END  K l 
END  J) 

GO  TO  STEP4I 

STEPS  I  PHASE40(ZIlJ*Sm-HRim/wm  t 
END  OF  PHASE4J 
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